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Abstract 

This  document  constitutes  the  final  project  report  for  Contract  #  F49620-99-C-0013 
titled  Tech  Sat  21  PHASED  ARRAY  SIGNAL  PROCESSING.  The  report  sununarizes 
the  theoretical  development  of  time-reversal  based  imaging  algorithms  for  locating 
and  tracking  ground  based  targets  from  multistatic  data  collected  from  unstructured 
phased  antenna  arrays  operating  above  the  ionosphere.  The  report  includes  a  number 
of  computer  simulated  examples  illustrating  the  use  of  the  algorithms  developed  in  the 
project. 
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1  Project  Summary 

Tech  Sat  21  has,  as  its  goal,  the  development  of  a  large,  sparse,  unstructured  3D  phased  array 
(Tech  Sat  21  array)  to  be  used  in  orbit  above  the  ionosphere  for  the  stuveillance  of  earth 
based  targets.  The  key  and  novel  ingredients  of  the  Tech  Sat  phased  array  system  is  that 
it  very  sparse  and  vmstructured;  i.e.,  the  phased  array  elements  are  not  rigidly  supported 
and  are  freely  orbiting  (uncontrollable)  in  three-dimensional  space.  These  requirements  of 
the  Tech  Sat  phased  array  system  require  that  new  and  novel  signal  processing  algorithms 
and  methods  be  developed  for  the  smrveillance  tasks  required  of  the  system.  The  research 
described  herein  is  directed  toward  developing  sudh  methods  and  algorithms  for  the  specific 
tasks  of  detecting  and  locating  (tracking)  moving  ground  targets  (MGT). 

The  approach  taken  in  the  project  is  based  on  the  well-known  concept  of  time-reversal 
imaging  [1,  2,  3]  which  will  be  explained  in  detail  in  the  body  of  the  report.  The  objectives 
of  the  research  were  to  develop  the  tmderlying  mathematical  and  physical  basis  for  (compu¬ 
tational)  time-reversal  imaging  for  the  Tech  Sat  project  and  to  develop  and  test  and  evaluate 
time-reversal  based  algorithms  for  MGT  detection  and  tracking.  Specific  areas  addressed  in 
the  research  included: 

•  An  in-depth  investigation  of  the  state-of-the-art  of  computational  time-reversal  imaging 
for  three-dimensional  (non-planar),  non-co-located  and  sparse  phased  arrays, 

•  Development  of  methods  to  incorporate  models  of  the  ionosphere  into  the  time-reversal 
based  algorithms. 
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•  Incorporation  of  state-of-the-art  signal  processing  schemes  such  as  MUSICS  into  the 
time-reversal  imaging  method, 

•  Development  of  methods  to  reduce  the  effect  of  clutter  on  time-reversal  imaging, 

The  major  accomplishments  of  the  project  were  in  the  development  of  a  MUSIC  algorithm 
for  MGT  location  and  the  generalization  of  the  usual  and  standard  methods  of  time-reversal 
imfl-ging  to  non-co-located  phased  array  systems  and  to  non-reciprocal  badcground  media 
such  as  the  ionosphere.  All  three  of  these  developments  have  or  are  in  the  process  of  being 
reported  in  the  literature  [4,  5,  6]  and  could  have  major  impact  on  the  Tedi  Sat  and  related 
projects.  In  addition  to  these  major  results  a  conceptual  scheme  for  the  complete  data 
processing  required  in  the  project  was  developed  and  a  number  of  computer  simulations 
were  completed.  These  and  other  results  from  the  project  are  reported  in  the  body  of  the 
report. 


2  Basic  Theory 


In  this  section  we  review  the  revamped  theory  of  computational  time  reversal  imaging  devel¬ 
oped  within  the  Tech-Sat  21  project.  One  of  the  principle  goals  of  the  Tedi-Sat  project  is  to 
employ  an  orbiting  set  of  radar  antennas  in  a  multistatic  mode  with  the  purpose  of  locating 
and  identifying  moving  ^oimd  targets  (MGT’s).  We  can  idealize  the  satellite  antenna  sys¬ 
tem  as  consisting  of  a  phased  array  system  that  interrogates  a  given  groimd  patch  containing 
an  MGT  as  illustrated^  in  figure  1.  In  our  analysis  we  will  assume  a  general  array  consisting 
of  Nr  receiving  antenna  elements  and  Nt  transmitting  elements,  all  of  which  interrogate  a 
given  groimd  patch  that  contains  a  total  of  M  targets.  Although  the  envisioned  Tech  Sat 
system  will  employ  co-located  transmit  and  receive  antenna  elements  the  generalize  theory 
is  necessary  to  treat  non-reciprocal  propagation  media  (such  as  the  ionosphere)  and  also  al¬ 
lows  the  results  obtained  in  the  project  to  be  employed  in  other  Air  Force  applications  (such 
as  imaging  of  ground  targets  from  phased  array  antennas  deployed  on  unmanned  aerial  ve¬ 
hicles  (UAV’s)).  The  effective  size  of  the  groimd  patch  is  determined  by  the  aperture  of  each 
antenna  element  (approximately  a  meter  diameter  in  Tech-Sat)  as  well  as  by  Doppler  and 
time-gate  processing.  The  Doppler  and  time-gate  processing  are  used  to  reduce  the  effect 
of  ground  clutter  and  are  required  pre-processing  steps  for  time-reversal  processing  to  be 
effective^.  We  illustrate  the  Doppler  preprocessing  steps  in  figure  2. 

The  pre-processed  data  (after  clutter  reduction)  at  frequency  co  forms  axx  Nr  xNt  matrix 
that  can  be  expressed  in  the  general  form: 


M  . 

m=i  ^ 


(1) 


^Standing  for  MUtiple-SIgnal  Classfication. 

^AU  figures  appear  at  the  end  of  the  report. 

^We  will  not  go  into  the  details  of  how  this  processing  is  done  since  it  is  standard  radar  theory  and  lies 
outside  of  the  main  goals  of  the'  research. 
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poterdich  chLact^izing  the  ground  targets  and  for  pemtrable  targets  are  given 

by  an  expression  of  the  form  i 

'  Om{r,u))  = 

where  P(r  w)  is  the  waven^tnber  of  the  target  and  fcg  the  wavenumber  of  the  background 
medium  We  emphasize  thait  the  above  model  makes  no  assumptions  about  the  9^^ 
of  the  phase  army  nor  of  reciprocity  of  the  background  medium.  It  is  also  exact  within  the 

framework  of  scalar  wave  scattering  theory. 

If  we  define  the  column  vectors 


Jfr(r,w)  lG(R;,r.a.),G(R5.r,w),  •  ■  •  ,G(RI,,.  r.«)f 

^(r,  lu)  =  (V»i  (r,  oi) ,  ^2(r,  w) ,  •  •  •  ,  (r,  <*^)  f 

we  can  rewrite  Eq.(l)  in  thes  symbolic  form 


(3a) 

(3b) 


/if(c^)  =  <i®r5r(r,w)Or„(r,w)V>^(r,u)  (4) 

in=l 

where  the  superscript  T  stajids  for  the  transpose  operation  and  where  we  have  used  the  r 
subscript  on  the  Green  function  vector  to  denote  that  it  corresponds  to  the  receiver  array 
The  data  matrix  K  =  is  called  the  multistatic  data  matrix  and  is  the  key  quantity 
that  is  employed  in  many  inverse  scattering  and  imaging  schemes.  In  Diffiraction  tomography 
(DT)  this  matrix  is  used  in  conjunction  with  a  slant  stacking  alpriihra  to  generate  plane  v^ve 
data  whidi  is  then  input  inti)  a  filtered  back  propagation  algorithm  tc>  generate  an  “image  of 
the  target  (the  object  profiles  O^).  Alternatively,  the  multistatic  matrix  can  be  employ 
to  perform  time-reversal  imaging  [1,  2]  whi^  is  especially  useful  for  locating  targets  t^t 
are  small  relative  to  a  wavelength,  even  if  the  targets  are  buried  o  mhomogeneous  n^ia 
whose  Green  functions  are  pot  precisely  known.  Finally,  we  mention  that  conventional  SAR 
imaging  does  not  use  the  eitire  multistatic  matrix  but  only  its  (Uagonal  eleinents  Ki,i  and 
is,  thus,  inherently  inferior  to  the  more  general  schemes,  such  as  t  ime-reversal  imaging,  that 
%mploy  the  entire  matrix. 


2.1  Born  Approxirtiation 

In  many  cases  of  interest  the  effects  of  the  targets  on  the  incident  wfive  field  ijjj  axe  negligible 
so  that  V'i  is,  to  a  good  approximation,  the  wave  field  that  is  ger  erated  by  the  j’th^tenna 
in  the  background  medium  without  any  targets  being  present.  This  is  called  the  Bom 
approximation^  and  is  employed  in  the  vast  majority  of  hnagiug  and  inversion  schemes. 

‘‘More  precisely,  this  is  called  the  distorted  Bom  approximaUon  since  .he  backpound  medium  is  not 
required  to  be  constant.  Howevet,  ‘we  wUl  interpret  the  Born  approximation  witain  this  more  gene^  context 
and  not  differentiate  between  uniform  (constant)  and  non-uniform  backgrounds  except  when  noted. 
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Another  eimpUacation  reeultein  cases  where  the  sensor  elements  (ant^»)  me  sm^  reMw 
t^e  lavetength.  In  this  cash  the  incident  waves  are  proportional  to  the  Green  fanctiona, 

I 


i.e. 


oj)  =  CiGir,  Rj-,  w) 

where  Ci  is  a  constant  (poesit>ly  frequency  dependent)  that  vre  will  take  to  unity  and  Rj 
rthe  trLsnutter  elemlt  location.  There  is  no  loss  of  generahty  ni  ^ng  this  ^-uP^n 
since  aU  sensor  array  constadts  can  be  easily  included  in  the  genei  al  formahsm.  Under  to 
assumption,  and  within  the  iom  approximation,  the  multistatic  resixinse  matnc  defined  m 

Elq.(4)  assumes  the  form 


M  - 


(5) 


where  Qr  is  the  Green  function  vector  associated  with  the  rece;  ver  array  and  defined  in 
Eq.{3a)and  •  ,G(r,Ri,„ai)f ,  (6) 

is  the  Green  function  vector  iassociated  with  the  transmitter  array.  Physically,  eaih  dement 
of  Pt  is  generated  by  a  transifiitter  antenna  point  and  propagates  fiom  that  transmitter  point 
into  the  target  re^on  while  is  generated  by  a  target  scattering  ]>oint  and  propagates  from 
that  target  scattering  point  Into  the  receiving  sensor  array.  This  is  illustrated  m  figure  3 

FVom  this  point  on  we  ^irill  employ  Eq.(5)  as  our  working  definition  of  the  multistatic 
response  matrix.  Thus,  we  ivill  assume  that  the  (distorted  wave)  Bom  approximation  can 
be  employed  and  that  the  Incident  wave  fields  generated  from  the  transmitter  array  are 
the  background  Green  functions.  This  second  assumption  is  done  only  for  convenience  and 
can  easily  be  removed.  For  ^the  sake  of  ease  of  notation  we  will  cJso  suppress  the  frequency 
variable  w  in  all  of  the  following  development  tvitii  the  understanding  that  all  equations  are 
in  the  frequency  domain  at  frequency  u. 

Finally,  we  emphasize  that  although  we  have  made  the  Bom  epproximation  we  have  not 
required  that  the  Green  fufictions  obey  a  reciprocity  condition.  For  example,  even  if  the 
transmitters  and  receivers  ate  co-  located  (as  is  presently  envisioned  for  Tech-Sat)  we  do  not 
require  that 

G(r,R5,a;)  =  G(R5,r,a;). 

Thus,  the  formulation  presented  here  can  be  employed  in  both  reciprocal  and  non-reciprocal 
badc^oimds  as  well  as  in  <4ses  where  the  transmit  and  receive  arrays  are  non-co-located. 


2.2  Time-reversal  Imaging 

In  the  general  case  where  all  of  the  transmitters  are  simultaneomly  activated  the  data  mea¬ 
sured  along  the  receiver  arrjay  are  given  by  an  expression  of  the  form 
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wh«-e  V  =  v(u})  is  the  lineax  array  of  output  voltages,  viewed  as  an  dimensional  column 
vector,  measured  at  the  receiver  array  terminals,  e  =  e(aj)  is  the  Nt  dimensional  column 
vector  of  apphed  excitations  to  the  set  of  transmitting  antennas,  and  K  is  the  Nr  x  Nt 
multistatic  matrix.  As  mentioned  earlier  we  will  not  explicitly  display  the  frequency  variable 
(j  in  subsequent  equations.  In  the  current  theory  of  time-reversal  imaging  the  object  profiles 
Om(r)  characterizing  the  targets  are  assumed  to  be  disjoint  profiles,  eadi  centered  at  a  spatial 
location  Xm  and  each  having  an  effective  size  that  is  small  relative  to  the  wavelength;  i.e., 

Om{r)  =  Om{r-Xm).  (7) 

The  goal  of  time-reversal  iTnaging  is  then  to  estimate  the  location  und  strength  of  each 
scatterer.  If  we  substitute  Eq.(7)  into  Eq.(5)  we  obtain 

d^rgr{T)Om{r  -  Xm)9f  (r) 

M 

«  (8) 

m=l 

where 

rm  =  J<^rOm{r)  (9) 

and  where  we  have  made  use  of  the  assumption  that  the  targets  are  small  relative  to  the 
wavelength.  The  quantities  thus  represent  effective  reflection  coefficients  for  the  targets 
and  the  goal  of  time-reversal  imaging  is  then  to  estimate  these  reflection  coeflScients  as  well 
as  the  target  locations  X^- 

2.3  Transmitter  and  Receiver  Array  Point  Spread  Functions 

Classical  coherent  imaging  from  arrays  is  performed  by  back  propagating  a  coherent  wave 
that  is  measured  across  the  array  into  the  space  from  which  it  propagated.  If  denotes 
such  a  wave  measured  along  the  receiver  array  then  the  back  propagation  of  ^  is  defined 
mathematically  via 

x(r)  =  f;Gr(r,R;WR;)  (10) 

where  x  is  the  back  propagated  wave  and  ^(Rp  the  measured  coherent  wave  across  the 
receiver  array.  The  back  propagated  wave  x(r)  is  the  “best”  estimate  of  the  measured 
coherent  wave  that  can  be  deduced  from  the  measured  data. 

Now  assume  that  the  incident  coherent  wave  to  the  receiver  is  the  Green  function 
G(R’‘,X)  resulting  from  a  source  location  at  X.  We  then  find  using  Eq.(lO)  that  the  back 


6 


propagated  field  corresponding  to  the  measured  Green  function  is 

Nr 

H,(r.X)  =  ^(r(r,R;)G{R;,X) 

=  5j(r)5r(X).  (11) 

The  back  propagated  field  Hr  is  a  function  both  of  the  field  point  at  which  it  is  evaluated 
as  well  as  the  source  point  location  X  and  is  the  best  image  of  the  source  point  that  can  he 
formed  from  measurement  of  the  Green  function  across  the  receiver  array.  This  quantity  is 
called  the  receiver  array  coherent  point  spread  function  (PSF).  In  a  similar  fashion  one  can 
define  the  transmitter  array  PSF  as 

H,(r,X)=  gl(r)g,(X).  (12) 

We  will  not  digress  into  the  importance  and  use  of  the  above  PSF’s  at  this  time  but  will 
shortly  encoimter  both  of  these  quantities  in  connection  with  time-reversal  imaging. 

2.4  SVD  of  the  multistatic  Matrix 

The  theory  of  computational  time-reversal  imaging  depends  on  the  ability  to  perform  a 
singular  value  decomposition  (SVD)  of  the  multistatic  data  matrix  K.  In  particular,  we 
consider  the  singular  system 

K:  Kej^ajVj  (13a) 

Xt :  ^  K^Vj  =  (Tjej  (13b) 

where  j  labels  the  singular  system  ej,Vj,{Xj.  The  normal  equations  for  this  system  are 

K^Kej  =  <T]ej,  (14a) 

KK^j  =  <t]vj.  (14b) 

The  singular  vectors  are  orthonormal  and  span  the  space  (7^*  while  the  singular 

vectors  are  orthonormal  and  span  the  space  There  are  a  total  of  min  {Nt,  Nr) 

singular  values  ctj  >  0. 

In  the  usual  theory  of  time-revers'al  imaging  [1,  2]  where  the  transmit  and  receive  arrays 
are  coincident  the  Nt  x  Nt  matrix 

T  =  K^K  (15) 

is  the  well-known  time-reversal  matrix.  It  is  seen  that  in  the  more  general  theory  developed 
here  that  thore  are,  in  fact,  two  time  reversal  matrices  Ti  =  K^K  as  well  as  T2  =  KK^.  Ti 
can  be  considered  to  be  a  “conventional”  time  reversal  matrix  for  a  single  coincident  sensor 
array  identical  to  the  transmit  array  and  T2  can  be  considered  to  be  a  “conventional”  time 
reversal  matrix  for  a  single  coincident  sensor  array  identical  to  the  receive  array. 
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If  we  substitute  the  expression  for  the  K  matrix  given  in  Eq.{8)  into  Eqs.(13a)  and  (16b) 
we  obtain 

M 

TtTi9rO^-rri)9t  == 

m=l 
M 

m—\ 

It  follows  from  the  above  equations  that  the  singular  vectors  Vj  having  non-zero  singular 
values  are  linear  combinations  of  the  receiver  array  Green  function  vectors  gr{X^)  while 
the  RiTigiilflT  vectors  ej  are  linear  combinations  of  the  complex  conjugates  of  the  transmitter 
array  Green  function  vectors  K  we  couple  these  observations  with  the  fact  that 

the  transmit  and  receiver  array  Green  function  vectors  are  linearly  independent  [4,  5]  we 
conclude  that  the  vector  space  spanned  by  the  transmit  (receive)  Green  function  vectors  is 
identical  to  the  vector  space  spanned  by  the  singular  vectors  e)  (vj)  having  non-zero  singular 
values  <Tj.  These  conclusions  form  the  basis  for  the  MUSIC  algorithm  discussed  below. 


(16a) 

(16b) 


2.5  Well-resolved  Targets 

An  important  special  case  occurs  when  the  two  sets  of  Green  function  vectors  are  orthogonal; 
i.e.,  when  the  following  two  equations  are  satisfied: 

pJ(Xm)^r(X„i')  =  |l5r(Xm)|P^m,m'  (17a) 

^|(XJpt(X„.0  =  IIPt(X^)||X.m'  (17b) 

where  Sm,m'  is  the  Kronecker  delta  function  and 

|k(X^)|P  =  5^(X^)^r(X^) 

\\gr{X^W  =  gl{Xm)gr{Xm) 

are  the  squared  norms  of  the  Green  function  vectors  evaluated  at  the  target  point  Xm-  If 
Eq.(17a)  holds  then  we  say  that  the  targets  are  well  resolved  by  the  receiver  array  while  if 
Eq.(17b)  holds  we  say  that  the  targets  are  well  resolved  by  the  transmitter  array.  When  both 
equations  hold  then  the  targets  are  well  resolved  with  respect  to  both  the  transmitter  and 
receiver  arrays. 

The  rational  for  the  above  terminology  is  apparent  if  we  simply  note  that  the  inner 
products  in  Eqs.(17a)  and  (17b)  are,  respectively,  the  PSF’s  Hr{Xm,  Xm’)  and  Ht{Xm,  Xm’). 
Thus,  for  example,  the  inner  product  gl{Xm)griXm')  is  the  image  of  a  point  target  located  at 
Xm'  formed  at  point  Xm  by  the  receiver  (receiver)  array.  An  entirely  analogous  interpretation 
can  be  given  the  inner  product  gt{Xm)gtiXm>)]  i-e-j  ss  the  image  of  a  point  target  located 
at  Xm’  formed  at  point  Xm  by  the  transmitter  array.  The  case  of  well  resolved  targets  thus 
corresponds  to  the  case  where  the  targets  are  sufficiently  well  separated  such  that  the  PSF  of 
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either  the  transmitter  or  receiver  array  does  not  significantly  overlap  any  target  other  than 
the  one  on  which  it  is  focused. 

For  the  class  of  targets  that  are  well-resolved  by  both  the  transmit  and  receive  array 
both  Eqs.(17a)  and  (17b)  hold.  In  this  case  it  is  easy  to  show  that  the  singular  system 
{ej,Vj,(Tj  >  0}  can  be  related  in  a  one-to-one  manner  with  the  M  <  min  {Nt,Nr)  isolated 
targets.  Indeed,  it  follows  at  once  firom  the  orthogonality  of  the  Green  function  vectors  that 
the  singular  vectors  having  non-zero  singular  values  for  well-resolved  targets  are  given  by 


„  9;(Xi) 

^  “  Il9.(x,)|| 

„  ■  »(Xi) 

’  ~  ll«r(X,)|| 


(18a) 

(18b) 

(18c) 


where  j  =  1,2,  Moreover,  the  non-zero  singular  values  rr,  are  given  by 

<rj  =  |ri|||s,(X,)||||s,(X,)||  (19) 

We  conclude  that  the  scatterer  strengths  are  computed  directly  from  the  singular  values  which, 
in  turn,  are  readily  computed  firom  the  measured  multistatic  data  matrix  K.  Moreover,  the 
singular  vectors  give  the  location  of  the  targets.  Indeed,  the  coherent  image  formed  using 
the  singular  vector  ej  will  generate  the  PSF  of  the  transmitter  array  centered  at  Xj  while 
the  coherent  image  formed  using  the  singular  vector  vj  will  generate  the  PSF  of  the  receiver 
array  centered  at  Xj. 

2.6  Non-well  resolved  TEU*gets  and  MUSIC 

One  of  the  primary  accomplishments  of  the  Tech-Sat  project  has  been  the  incorporation 
of  MUSIC  into  computational  time-reversal.  MUSIC  allows  non-well  resolved  targets  to  be 
detected  and  located.  The  general  idea  is  that  although  the  set  of  transmit  and  receive 
Green  function  vectors  are  not  orthogonal  in  the  general  non-well  resolved  case  they  are  still 
linearly  independent  (This  result  is  proven  in  [4])  and,  as  indicated  earlier,  this  impUes  that 
the  vector  space  spanned  by  the  transmit  (receive)  Green  function  vectors  is  identical  to  the 
vector  space  spanned  by  the  singular  vectors  Cj  (vj)  having  non- zero  singular  values  Cj.  It 
then  follows  that 


5;^  ets,-(X„.)  =  0  (20a) 

<7^=0 

^uV(X„.)=0  (20b) 

where  the  sums  are  over  all  singular  vectors  having  zero  singular  value  and  where  X^  is  any 
target  location.  It  is  important  to  note  that  the  above  equations  hold  independent  of  whether 
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the  targets  are  well  resolved  or  not.  The  orthogonality  conditions  given  in  Eqs.(20)  is  the 
underlying  key  to  MUSIC.  In  particular,  we  introduce  the  steering  vectors 

^t(Xp,a;)  =  [G(Xp,Rl,w),G(Xp,R4,a;),---  ,G(Xp,R5^,,a;)p,  (21a) 

gr{Xp,u)  =  [G(Rl,Xp,a;),G(R5,Xp,a;),---  ,G(R^^,Xp,a;)f  (21b) 

where  Xp  is  a  test  target  location  vector  arid  can  assiune  any  value  within  the  region  occupied 
by  the  targets.  It  then  follows  that  when  the  test  target  location  vector  coincides  with  an 
actual  target  location;  i.e.,  when  Xp  =  Xm  then 

<7j=0 

Y,  "iSrCXp)  =  0. 

Orj=0 


MUSIC  uses  the  pseudo-spectrum 


P(Xp)  = 


E,,M,H9f(Xp)  +  <'V{X,)lP 


to  determine  the  target  locations.  It  follows  from  the  analysis  presented  above  that  V{Xj,) 
will  become  large  (infinity  in  the  ideal  case)  at  all  target  locations  (Xp  =  X^)  aJid  will 
be  sTna.11  otherwise.  A  number  of  computer  simulated  examples  of  the  use  of  the  pseudo¬ 
spectrum  in  the  special  case  where  the  transmit  and  receive  arrays  are  co-located  are  pre¬ 
sented  in  reference  [4]  .  The  case  of  non-co-located  arrays  are  treated  in  reference  [5].  We 
present  simulations  for  the  case  of  co-located  arrays  in  the  following  section. 


3  Near  Field  Computer  Simulations  for  Tech  Sat  21 

We  i>erformed  a  number  of  computer  simulations  to  test  and  evaluate  the  computational 
time  reversal  algorithms  developed  dining  the  course  of  the  project  and  outlined  above.  The 
simplest  of  these  simulations  were  performed  for  the  case  of  near  field  imaging  where  the 
objects  were  located  in  the  near  vicinity  of  the  phased  array.  Although  this  is  not  directly 
relevant  to  Tech  Sat  it  is  relevant  to  a  number  of  Air  Force  and  DOD  applications  such 
as  ground  imaging  from  Unmanned  Aerial  Vehicles  (UAV’s).  These  initial  simulations  also 
serve  as  a  test  of  the  basic  performance  of  the  algorithms  in  an  idealized  setting.  To  avoid 
discontinuities  in  the  presentation  and  to  avoid  confusion  we  have  included  all  simulation 
figures  at  the  very  end  of  the  document. 

The  simulations  were  coded  using  MATLAB  and  assume  a  two-dimensional  geometry 
where  the  sensor  elements  are  line  sources  and  the  targets  are  line  targets  all  embedded  in 
a  homogeneous  background  and  perpendicular  to  a  single  plane  and  all  multiple  scattering 
is  ignored.  FVom  a  mathematical  point  of  view  we  can  consider  the  targets  iand  sensor 
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elements  as  points  in  the  two-dimensional  space  of  the  plane;  i.e.,  as  point  locations  in  F^. 
In  the  simulations  we  have  hmited  our  attention  to  the  case  of  co-located  arrays  which  are 
currently  envisioned  for  Tech  Sat.  In  these  cases  there  is  only  one  Green  function  vector 
which  we  can  take  to  be  the  transmitter  Green  function  vector  defined  in  Eq.(3a).  In  the 
homogeneous  backgroimd  case  the  appropriate  Green  function  G(r,r')  is(  given  by 

G(r,r')  =  ^fro(fc|r-r'|)  (23) 

where  Hq  is  the  zero  order  Hankel  function  of  the  first  kind.  The  Green  function  vectors  are 
then  given  by 

.  9m  =  gr{^m)  =  {Ho{k\Ri  -  X^|)} 

=  [Ho{k\Ri  -  X^l),  Ho{k\R2  -  X^l), . . . ,  fl-o(fe|Rjv  -  X^|)F,  (24) 

where  we  have  dropped  the  unessential  constant  —i/4  and  where  we  have  denoted  the  sensor 
element  locations  of  the  phased  array  by  Rj,  /  =  1, 2,  •  •  •  ,  TV.  Note  that  we  use  the  shorthand 
notation  gm  to  denote  the  array  Green  fimction  vectors  gr  evaluated  at  the  scatterer  locations 
X^;  i.e.,  gm  =  grp^m)-  In  terms  of  the  pr  the  multistatic  response  matrix  (MSR)  K  is  given 
by 

K  =  '^rmgm9m-  (25) 

m 

In  the  simulations  we  employed  a  basic  image  space  which  is  a  rectangular  grid  repre¬ 
senting  a  section  of  the  x,  z  plane  with  z  being  depth  (pointing  down)  and  x  lateral  location 
(with  positive  x  directed  to  the  right).  All  dimensions  are  relative  to  the  wavelength  which 
is  taken  to  be  unity.  The  simulations  used  an  image  grid  spacing  parameter  5x  —  Sz  equal 
to  a  quarter  wavelength  Sx  =  A/4  and  an  uniformly  spaced  linear  array  of  firom  five  to 
nine  elements  (depending  on  the  simulation)  located  along  the  line  z  =  0  and  with  adjacent 
sensor  element  spacing  var3dng  from  a  half-wavelength  to  up  to  sixteen  wavelengths,  again 
depending  on  the  simulation.  We  also  used  from  two  to  four  targets  located  at  the  same 
z  (depth)  coordinate  (sixteen  wavelengths)  but  with  variable  spacing  along  the  x  direction. 
In  some  of  the  simulations  we  added  additive  noise  to  the  computed  multistatic  response 
matrix  (MSR  matrix)  K  according  to  the  model 

K  =  K  +  (26) 

where  K  is  the  noisy  MSR  matrix,  .4  is  an  independent  Gaussian  variable  with  variance 
proportional  to  the  peak  amphtude  of  the  ideal  (computed)  MSR  matrix  K  and  W  is  an 
independent  Gaussian  variable  with  unit  variance.  Thus,  the  noise  model  includes  both 
amphtude  and  phase  fluctuations.  All  parameters  used  in  any  given  example  are  listed  in 
the  figure  captions. 
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3.1  Computing  the  multistatic  Response  Matrix 

The  multistatic  response  matrix  K  is  computed  using  Eki.(25)  with  the  Qm  equal  to  the 
receiver  Green  function  vectors  gr  evaluated  at  r  =  X,„.  As  an  initial  example  we  computed 
the  MSR  matrix  for  the  case  of  two  targets  (M  =  2)  both  located  at  a  depth  of  sixteen 
wavelengths  for  a  linear  sensor  array  of  nine  elements  that  was  centered  over  the  image  grid. 
In  this  first  simulation  we  used  a  half-wavelength  spacing  for  the  sensor  elements,  no  additive 
noise,  and  a  target  separation  of  eight  wavelengths  along  the  x  direction,  with  one  target 
located  exactly  at  the  mid-point  of  the  array  and  the  other  displaced  by  eight  wavelengths 
in  the  positive  x  direction.  The  scattering  strengths  of  the  two  targets  were  equal  (to 
unity).  The  Coherent  Point  Spread  Function  (CPSF)  of  a  regularly  spaced  sensor  array  of 
length  a  will  have  an  effective  lateral  width  S  at  distance  2:  of  approximately 

S  =  zX/a 

which  using  z  =  16A,  a  =  4A  yields  S  =  4A.  This  is  equal  to  half  of  the  spacing  between  the 
two  targets  so  this  example  corresponds  to  two  well-resolved  targets. 


3.2  Computing  the  Eigenvalues  and  Eigenvectors  of  the  time- 
reversal  matrix 

The  time-reversal  matrix  T  is  computed  directly  from  the  MSR  matrbc  according  to  the 
equation 

T  =  K^K  =  r^rm>gl^\n9m'9Z.> 

mm* 

=  (27) 

m  m' 


where 


Afn,TO'  —  ’^mTfn'9m9tn'‘ 


We  used  a  standard  eigenvector/eigenvalue  solver  in  MATLAB  to  compute  the  eigenvalues 
and  eigenvectors  of  the  time-reversal  matrix  T  in  the  simulation.  Since  there  are  two  targets 
there  are  two  non-zero  eigenvalues  with  associated  signal  space  eigenvectors.  We  show  in 
Fig.  5  a  plot  of  the  eigenvalues  and  unwrapped  phase  of  the  eigenvectors  computed  from 
the  MSR  matrix  and  for  the  simulation  parameters  given  above.  We  see  for  this  case  that 
there  are  two  non-zero  eigenvalues  corresponding  to  the  two  well  resolved  targets  that  are 
separated  by  eight  wavelengths.  We  also  show  in  the  bottom  figure  the  phase  of  the  complex 
conjugate  Green  function  vectors  g’^  which  are  seen  to  be  nearly  identical  (to  within  an 
additive  integral  multiple  of  27r)  to  the  phase  of  the  two  eigenvectors.  This  is  in  agreement 
with  the  fact  that  the  complex  conjugate  of  the  Green  function  vectors  are  proportional  to 
the  eigenvectoirs  of  the  time-reversal  matrix  in  the  case  of  well-resolved  targets  [4].  We  show 
only  the  phase  of  these  quantities  since  the  phase  of  the  field  is  much  more  important  than 
amplitude  (intensity)  in  imaging  and  target  location  estimation. 
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3.3  Forming  the  conventional  time-reversal  image 

The  time-reversal  image  is  formed  using  Eq.(lO).  For  the  well  resolved  target  case  the  image 
field  V’m(r)  is  proportional  to  the  Coherent  Point  Spread  Function  (CPSF)  of  the  sensor 
array.  This  is  well  approximated  in  the  first  simulation  where  the  targets  are  separated 
by  8A.  The  magnitude  of  the  focused  fields  from  the  two  eigenvectors  for  this  case  (whose 
phases  are  shown  in  Fig.  5)  are  shown  in  Fig.  6.  We  have  superimposed  an  X  at  the  actual 
target  locations  to  show  that  the  image  fields  tend  to  peak  in  the  vicinity  of  the  targets. 
We  have  also  shown  for  comparison  the  magnitude  of  the  fields  generated  from  the  complex 
conjugate  Green  function  vectors  which  are  the  eigenvectors  in  the  ideal  case  where  the 
targets  are  perfectly  resolved  and  whose  phase  distributions  along  the  sensor  elements  are 
also  shown  in  Fig.  5.  It  is  clear  firom  the  figure  that  the  image  intensities  for  both  the  imaged 
eigenvectors  and  the  imaged  Green  function  vectors  are  very  close  as  can  also  be  inferred 
from  the  fact  that  their  phase  distributions  along  the  sensor  elements  are  almost  identical 
as  shown  in  Fig.  5. 

It  is  clear  firom  Fig.  6  that  the  images  generated  firom  the  eigenvectors  of  the  time- 
reversal  matrix  provide  a  good  indicator  of  the  target  locations.  The  reasons  for  this  are  that 
the  targets  are  well-revived  so  that  their  separate  image  fields  do  not  overlap  significantly. 
However,  we  will  see  later  examples  where  the  targets  are  not  well-resolved  so  that  the  images 
generated  by  the  eigenvectors  of  the  time-reversal  matrix  do  not  provide  good  indication  of 
the  target  locations.  Moreover,  even  in  this  example  where  the  targets  are  well  separated  the 
depth  resolution  generated  by  the  classical  time-reversal  images  is  not  very  good.  Indeed, 
a  close  inspection  of  the  images  in  Fig.  3  shows  that  the  maximums  do  not  occur  at  the 
target  locations  but  rather  in  the  immediate  vicinity  of  the  array  (see  gray  scale  bar).  The 
reason  for  this  is  that,  as  was  mentioned  earlier,  the  longitudinal  resolution  associated  with 
the  CPSF  is  much  less  than  is  the  horizontal  resolution:  a  fact  that  is  clear  firom  Fig.  6. 


3.4  Computing  the  Pseudo-spectrum 

The  MUSIC  algorithm  was  implemented  using  Elq.(22)  with  the  denominator  computed  using 
two  different  methods  that  are  described  in  [4].  Note  that  in  this  case  where  the  transmitters 
and  receivers  are  co-located  that  the  pseudo-spectrum  reduces  to  the  simplified  form 


r>(x,)  = 


^<r,=0 


\e]9p? 


where  Qp  =  pm|m=p  =  Pr(Xp)  with  Xp  equal  to  the  location  of  a  test  target.  In  all  of 
the  simulations  the  two  methods  of  computing  the  denominator  gave  essentially  equivalent 
results  so  we  only  show  the  result  obtained  by  the  direct  method  of  projecting  gp  onto  the 
noise  subspace.  In  the  first  simulation  discussed  so  far  the  targets  are  well-resolved  so  that 
each  eigenvector  is  associated  with  one  of  the  two  targets  and  the  conventional  method  of 
imaging  the  eigenvectors  illustrated  in  Fig.  6  gives  good  indication  of  target  location,  at  least 
as  regards  the  lateral  (x)  location.  The  computed  pseudo-spectrum  for  this  example  is  shown 
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in  Fig.  7  while  Fig.  8  shows  the  pseudo-spectrum  superposed  on  top  of  the  conventional  time- 
reversal  images.  Also  show  in  text  are  the  estimates  of  target  location  obtained  by  simply 
finHing  the  TnaxininTn  of  the  pseudo-spectrum.  As  expected  in  the  absence  of  noise  perfect 
target  location  estimation  is  obtained. 

3.5  Further  Examples 

As  a  second  example  we  repeated  the  first  simulation  but  added  noise  equal  to  20%  of 
the  peak  value  of  the  magnitude  of  the  (noise  free)  MSR  matrix  and  also  increased  the 
spacing  between  adjacent  sensor  elements  to  two  wavelengths.  In  Fig.  9  we  show  the  plot 
of  the  eigenvalues  as  well  as  the  phase  of  the  eigenvectors  and  complex  conjugate  Green 
functions.  Note  that  the  noise  has  the  effect  of  adding  non-zero  eigenvalues  beyond  the 
first  two  dominant  ones.  Also  note  that  the  phases  of  the  conjugate  Green  functions  are 
no  longer  equal  to  the  phases  of  the  first  two  eigenvectors  due  to  the  effect  of  the  increased 
sensor  separation  and  also  the  additive  noise.  In  Fig.  10  we  show  the  images  generated  from 
the  two  dominant  eigenvectors  together  with  those  generated  using  the  complex  conjugate 
Green  function  vectors.  It  is  seen  that  these  images  possess  a  complicated  spatial  structme 
with  multiple  lobes  making  it  difficult  to  located  the  two  targets.  In  Fig.  11  we  show  the 
pseudo-spectrmn  which,  although  noisy,  returns  exact  estimates  of  the  target  locations. 


4  Tech  Sat  Simulations 

The  Tech-Sat  simulations  included  the  presence  of  the  ionosphere  and  are  also  performed 
in  the  far  field  of  the  targets.  Both  of  these  aspects  of  the  problem  complicate  the  imaging 
problem  and  the  simulations  of  the  time-reversal  based  imaging  algorithms  discussed  in  the 
report.  As  in  the  near  field  simulations  presented  in  the  previous  section  we  consider  the 
simplified  case  of  two  space  dimensions  where  the  horizontal  x  axis  is  parallel  to  the  locally 
flat  earth  smfface  and  the  vertical  z  axis  has  its  origin  at  a  height  h  above  the  earth  with 
positive  z  pointing  downward.  A  phased  antenna  array  consisting  of  N  antenna  elements 
each  with  a  uniform  apertmre  of  diameter  d  is  assumed  to  be  deployed  along  the  x  axis  at 
2  =  0  (at  a  height  h  above  the  earth  surface).  Intervening  between  the  phased  anterma  array 
and  the  earth  is  the  ionosphere  which  is  modeled  as  a  thin  phase  screen  located  at  altitude 
I  above  the  earth.  The  thin  phase  screen  model  is  illustrated  in  fig.  4.  The  objective  of  the 
simulation  is  to  illustrate  the  use  of  the  time-reversal  based  processing  schemes  outlined  in 
the  previous  sectibns  of  the  report  for  locating  one  or  more  moving  ground  targets  (MGT) 
from  the  multistatic  data  matrix  acquired  at  a  single  frequency  by  the  orbiting  array.  We  will 
assume  that  all  data  collection  is  performed  over  a  short  period  of  time  in  which  everything 
is  “frozen”  so  that  variations  in  anterma  location,  target  location,  ionosphere  parameters, 
etc.  can  be  ignored.  The  assumption  of  a  frozen  backgroimd  and  data  acquisition  sjrstem  is 
reasonable  within  the  context  of  Tech  Sat  due  to  the  planned  use  of  frequency  multiplexing 
in  the  acquisition  of  the  multistatic  data  matrix.  Thus,  in  particular,  this  quantity  will  be 
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obtained  by  simultaneously  exciting  orthogonal  signals  at  all  the  transmitter  locations  and 
not  by  sequentially  exciting  the  transmitter  elements  in  time. 

We  will  take  the  nominal  altitude  h  of  the  Tech  Sat  array  to  be  400  kilometers  (km)  and 
place  the  thin  phase  screen  model  of  the  ionosphere  at  1  =  /i/2  =  200  km.  The  wavelength 
of  the  radiation  is  talcen  to  be  1  cm— 10“^  m  and  we  will  assume  that  the  array  consists 
of  i\r  =  5  antennas  aligned  along  the  x  axis  at  z  =  0  (height  h)  and  each  of  which  has  an 
effective  aperture  diameter  of  one  meter.  Due  to  the  large  propagation  distances  and  short 
wavelength  we  can  employ  the  following  FVaunhofer  approximation  for  the  field  radiated 
above  the  ionosphere  by  an  antenna  element  located  at  o:*  =  (a:  =  Xk,  2  =  0): 

$(r,  otk)  =  y  ^=^e*^^*"**')*sinc  \—{x  -  rrjt)] ,  (28a) 

where  A^)  =  27r/ A  is  the  free  space  wavenumber  and 

,  .  sinnx 
smc  (x)  = - 

TTX 

is  the  sine  function.  The  field  O  is  radiated  from  each  antenna  element  down  to  the  thin 
phase  screen  where  it  is  spatially  modulated  by  a  transmission  function  of  the  general  form 

r(a;)  = 

where  the  amplitude  Aq  and  period  Lo  are  parameters  of  the  thin  phase  screen  model.  In 
the  simulations  we  took  .4o  =  1  and  Lq  =  30  m.  More  general  models  for  the  phase  screen 
can  be  employed  and  are  easily  incorporated  into  the  Matlab  a)de  which  is  available  from 
the  author’s  web  site. 

The  thin  phase  screen  transmission  function  T{x)  multiplies  the  antenna  field  ^(r,  a*) 
at  2  =  Z  =  h/2=200  m  to  generate  the  boundary  v^ue  field  x{r,otk)  =  T{x)^{r,ak)  which 
is  then  propagated  to  ground  using  a  discretized  Kirdioff  approximation 


$(r,o*)  =<5x  ®^>Wc  [3^(2^ -  ^i)] 


(28b) 


where  oik)  are  the  sample  values  of  the  boundary  value  field  at  the  sample  points 
=  (jSXf  1),  Az  —  z—l  is  the  propagation  distance  from  the  phase  screen  and  the  maximum 
index  J  is  selected  large  enough  to  encompass  the  major  lobe  of  the  radiation  pattern  from 
each  anteima  element.  Eq.(28b)  can  be  regarded  to  be  a  superposition  of  the  fields  generated 
by  a  linear  array  of  adjacent  small  antenna  elements  having  diameters  Sx  and  boundary  value 
fields  x(rj » <**)•  In  the  simulations  we  took  Sx  =  d=l  m  since  the  distance  of  the  phase  screen 
from  the  phased  array  was  sufficiently  large  that  this  value  of  Sx  was  well  below  the  required 
Nyquist  sampling  interval. 

$(r,ajk)  is  the  total  field  radiated  by  an  antenna  centered  at  a*  in  the  presence  of  the 
thin  phase  screen®.  We  will  also  have  need  of  the  wavefield  which  is  the  output 

®We  assume  that  reflections  by  the  ionosphere  can  be  ignored  or  time  gated  out  by  the  Tech  Sat  antenna 
elements  so  that  the  reflected  wave  component  of  $  is  effectively  zero. 
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from  the  k  antenna  element  due  to  a  point  somrce  located  on  the  ground  at  x^-  This  field 
is  computed  using  the  model  Eq.(28a)  for  the  up  going  wave  fi:om  a  source  point  at  and 
then  the  Kirchoff  model  Eq.(28b)  for  propagation  up  to  the  antenna  element  located  at  a*;. 
However,  due  to  reciprocity,  it  is  easily  verified  that® 

^{ock,  Xm)  =  ^(Xm,  OCk).  (29) 

The  field  ^(r,Ofc)  plays  the  role  of  the  Green  function  G(r,aA:)-'  In  place  of  the  Green 
function  vector  we  now  have  the  “antenna  vector” 

<^(r)  =  [$(r,  CKi),  $(r,  0:2),  •  •  •  ,  ^(r,  ojat)]^,  (30a) 

where  ^  is  defined  in  Eqs.(28).  The  multistatic  data  matrix  is  given  by 

M 

/iC  =  5^Tm<^(Xm)0^(Xm),  (30b) 

m=l 

which  follows  at  once  firom  the  reciprocity  condition  Eq.(29). 

In  some  of  the  simulations  uncorrelated  noise  was  added  to  the  output  firom  the  antenna 
elements  according  to  the  model 

kj^k  =  Kj,k  +  Ai  exp  iXj,k 

where  K  is  the  noisy  matrix  and  Ai  is  a  real  constant  which  we  took  to  be  either  zero 
or  0.1  times  the  maximum  value  of  the  noise  firee  K  matrbc.  The  quantity  Xj^k  is  a  zero 
moan,  uncoixelated  anH  uniformly  distributed  random  matrix  over  [— tt,  tt].  We  note  that 
the  assumption  of  point  targets  is  very  reasonable  in  the  current  application  in  view  of  the 
very  large  propagation  distances  involved  even  though  the  wavelength  is  quite  small.  In 
particular,  the  incident  field  ^{r,ak)  from  each  antenna  element  evaluated  on  the  groimd 
will  be  effectively  plane  over  any  reasonably  sized  ground  vehicle  with  the  result  that  the 
point  scattering  model  imderlying  the  analysis  presented  in  the  paper  should  be  valid. 

4.1  Simulation  Results 

The  Matlab  code  apOl.m  essentially  tests  the  forward  simulation  model  Eqs.(28).  Among 
other  computations,  this  code  computes  the  field  intensity  over  the  earth  surface  (z=400 
km)  with  and  without  the  presence  of  the  phase  screen.  An  example  is  presented  in  fig.  12 
which  shows  the  field  intensity  at  the  top  of  the  phase  screen  as  well  as  over  the  ground 
with  and  without  the  phase  screen  being  present.  The  parameter  values  for  this  simulation, 
and  all  other  simulation  examples,  are  given  in  Table  1.  In  this  example  we  have  used  a 
uniform  linear  array  oiN  =  h  elements  equally  spaced  over  a  half  kilometer  and  all  of  which 

®In  all  the  Matlab  codes  discussed  later  and  used  in  the  simulations  the  up  going  wave  field  ®  is  actually 
computed  using  E)qs.(28)  and  the  reciprocity  condition  can,  in  fact,  be  verified  using  these  codes. 
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were  simultaneously  excited  by  unit  amplitude,  in-phase  pulses  to  simulate  a  plane  wave 
propagating  parallel  to  the  z  axis.  The  middle  and  top  plots  in  the  figure  clearly  show  the 
grating  lobes  that  are  generated  for  this  case.  The  grating  lobe  periods  are  given  in  the  table 
where  the  first  number  is  for  the  top  of  the  phase  screen  and  the  second  for  the  earth  Surface. 
It  is  also  apparent  firom  the  bottom  plot  that  the  grating  lobes  are  randomized  due  to  the 
presence  of  the  phase  screen.  The  code  apOl.m  can  also  be  used  to  test  the  accuracy  of  the 
simple  Kirchoff  model  Eq.(28b)  for  generating  the  incident  field  below  the  phase  screen.  In 
particular,  the  reader  can  verify  using  the  code  that  the  field  amplitude  on  the  earth  surface 
obtained  using  the  direct  model  Eq.(28a)  with  z  =  400  km.  agrees  closely  with  the  field 
amplitude  computed  using  the  composite  model  Eq.{28b)  with  the  phase  screen  amplitude 
Aq  set  equal  to  zero;  i.e.,  with  T{x)  =  1. 

The  code  apOB.m  allows  the  user  to  use  up  to  five  targets  distributed  over  a  half  kilometer 
on  the  earth  surface  and  computes  the  multistatic  data  matrix  K  and  the  eigenvectors  and 
eigenvalues  of  the  time-reversal  matrix  T  =  K*K.  These  eigenvectors  are  equal  to  the 
singular  vectors  vj  =  Uj  that  are  used  in  the  computation  of  the  classical  time-reversal 
image  field  via  the  process  of  back  propagation  as  discussed  earlier  in  the  report.  For  the 
case  of  a  single  target  or  a  set  of  “well  resolved  targets”  the  time-reversed  image  field  is  the 
PSF  of  the  antenna  array  centered  at  the  target(s)  locations.  As  discussed  earlier  the  PSF 
yields  the  “best”  image  of  the  target  location  that  can  be  generated  firom  the  observed  data 
(the  K  matrix). 

We  show  in  fig.  13  plots  of  the  intensity  of  the  time-reversed  field  over  the  groimd  for  a 
single  target  located  at  the  center  =  0  of  a  half  kilometer  of  image  space  using  the  same 
set  of  parameters  employed  in  the  first  example.  The  figure  shows  the  results  for  three  cases: 
(i)  when  no  phase  screen  is  present,  (ii)  when  a  phase  screen  is  present  and  (iii)  when  a  phase 
screen  is  present  but  the  image  field  is  generated  by  using  the  free  space  model  Eq.(28a). 
All  three  cases  were  computed  by  first  generating  the  K  matrix  and  then  back  propagating 
the  eigenvector  corresponding  to  the  largest  eigenvalue.  The  third  case  corresponds  to  using 
the  free  space  anteima  vector  generated  using  the  free  space  propagation  model  Eq.(28a) 
with  z  =  400  km  in  the  computation  of  the  image  field  and  is  of  interest  due  to  the  fact 
that  an  adequate  model  of  the  thin  phase  screen  (the  backgroimd  Green  function)  may  not 
be  available.  Thus,  although  the  exact  eigenvector  of  the  K  matrix  can  be  computed  the 
step  of  image  formation  via  back  propagation  may  not  be  exactly  implementable  due  to 
imperfect  knowledge  of  the  background  Green  function.  As  indicated  in  the  figure  the  use 
of  free  space  back  propagation  shifts  the  location  of  the  target  but  does,  nevertheless,  give  a 
rough  estimate  of  its  position. 

We  show  in  fig.  14  results  for  the  case  of  an  irregular  (unstructured)  antenna  array  having 
center  locations  as  indicated  in  Table  1.  We  have  again  used  a  single  target  located  at  the 
center  ajm  =  0  and  again  show  the  three  cases  of  no  phase  screen,  phase  screen,  and  phase 
screen  data  but  free  space  back  propagation.  Note  the  elongated  period  in  the  grating  lobe 
structure  compared  with  that  displayed  in  fig  13  due  to  the  smaller  minimum  separation  of 
antenna  elements  of  the  phased  array. 

It  is  apparent  from  figures  13  and  14  that  the  conventional  time  reversed  image  field  yields 
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good  estimates  of  the  target  location  in  the  absence  of  noise  and  clutter  if  the  background 
Green  function  (or  antenna  vector)  is  known  and  fair  estimates  if  this  quantity  is  not  known 
and  free  space  back  propagation  is  employed.  In  fig  15  we  show  the  effect  of  multiple  targets 
(clutter)  on  the  image  formation  process.  This  example  was  generated  with  a  phase  screen 
present  using  the  code  apOS.m  and  the  parameters  given  in  Table  1.  The  example  used  three 
targets  having  equal  scattering  strengths  =  1  located  at  Xm  =  -200,0,  and  -f  100  meters 
and  illustrates  the  difficulty  of  associating  individual  targets  from  the  raw  back  propagated 
images  even  in  the  absence  of  additive  noise.  The  problem  is  exacerbated  in  the  absence  of 
the  phase  screen  where  the  grating  lobes  make  target  location  even  more  difficult. 

4.2  MUSIC 

The  Matlab  code  apO^.m  implements  the  MUSIC  algorithm.  We  applied  this  code  to  the 
same  target  and  antenna  geometry  employed  in  the  example  given  in  fig.  15  and  present  the 
results  in  figs.  16  and  17,  The  results  shown  in  fig.  16  correspond  to  no  phase  screen  with  the 
top  plot  being  the  pseudo-spectrum  for  the  no  noise  case  and  the  bottom  for  an  additive  noise 
having  an  ampUtude  coefficient  Ai  =  .1  corresponding  to  a  signal  to  noise  ratio  maxJRr/maxW 
of  ten  where  W  is  the  additive  noise  term.  This  same  noise  amplitude  and  signal  to  noise 
ratio  was  employed  in  all  of  the  simulations  having  noise  present.  As  mentioned  in  our 
discussion  on  grating  lobes  the  multistatic  data  matrix  is  not  exactly  periodic  due  to  its 
dependence  on  the  square  of  the  target’s  location  so  that  in  the  absence  of  noise  it  should 
be  possible  to  perfectly  locate  a  target  even  in  the  presence  of  these  lobes.  Although  this  is 
not  easily  accomplished  by  direct  viewing  of  the  classical  time-reversed  (back  propagated) 
image  fields  the  target  locations  are  immediately  and  automatically  generated  by  the  pseudo¬ 
spectrum  as  is  apparent  from  the  top  plot  in  fig.  16,  However,  when  even  a  small  amount 
of  noise  is  added  a  nmnber  of  false  peaks  in  the  pseudo-spectrum  appear  that  are  due  both 
to  additional  clutter  like  singular  vectors  in  the  noise  subspace  H  as  well  as  to  the  basic 
ill-posedness  of  the  inverse  problem  in  the  presence  of  the  grating  lobes.  This  is  illustrated 
in  the  bottom  of  the  figure  where  a  small  amount  of  additive  noise  with  the  noise  amplitude 
coefficient  Ai  =  .1  has  been  added. 

In  fig.  17  we  show  the  results  obtained  for  the  same  set  of  simulation  parameters  used  for 
fig.  16  but  with  a  phase  screen  present.  The  randomizing  effect  of  the  presence  of  the  phase 
screen  on  the  grating  lobes  is  apparent  from  the  bottom  plot  which,  on  comparison,  with  the 
bottom  plot  of  fig.  17  is  seen  to  yield  much  more  robust  estimates  of  the  target  locations.  We 
performed  the  same  simulation  but  used  the  free  space  steering  vector  and  show  the  results 
in  fig.  18.  It  is  seen  that  the  free  space  steering  vector  does  not  give  accurate  estimates  of 
target  location  but  does,  nevertheless  yield  estimates  that  are  in  the  general  ball  park  of  the 
correct  locations. 

The  final  example  presented  in  fig,  19  is  of  two  close  targets  in  the  presence  of  clutter 
and  additive  noise.  This  example  is  important  since  actual  field  data  will  contain  a  nmnber 
of  clutter  targets  as  well  as  noise.  Both  examples  assumed  the  presence  of  a  phase  screen 
with  the  top  plot  corresponding  to  the  no  noise  case  while  the  bottom  plot  included  an 
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additive  noise  with  noise  amplitude  Ai  =  .1.  The  two  plots  clearly  show  the  robustness  of 
the  MUSIC  algorithm  and  also  indicate  its  ability  to  “super-resolve”  close  targets  that  would 
not  be  resolvable  using  standard  field  back  propagation.  Indeed,  the  two  dominant  targets 
in  this  example  are  separated  by  five  meters  which  is  roughly  equal  to  the  classical  resolution 
A  =  A/i/.S  =  4  m  of  a  densely  packed,  uniform  antenna  array  one-half  kilometer  in  length. 

5  Summary 

This  project  had  as  its  primary  goal  the  development  of  theory  and  associated  algorithms  for 
locating  gromid  based  targets  using  the  proposed  Tech  sat  phased  array  anteima  system.  To 
this  end  a  theory  based  on  classical  time-reversal  imaging  [1,  2]  was  developed  and  coded  in 
Matlab  based  algorithms.  The  basic  theory  has  been  submitted  for  publication  [4,  5,  6]  and 
is  summarized  in  this  final  project  report.  The  developed  theory  is  apphcable  to  Tech  Sat  as 
weU  as  to  general  problems  that  require  imaging  of  targets  obscured  by  an  inhomogeneous 
background  medimn  such  as  the  ionosphere  or  fohage. 

The  theory  and  algorithms  developed  in  the  project  assume  measurement  of  the  multi¬ 
static  data  matrbc  K  =  {Kj,k}  at  a  single  temporal  frequency  w  by  a  phased  array  antenna 
and  uses  the  S  VD  of  this  matrbc  as  well  as  the  Green  function  of  the  backgrovmd  medimn  to 
generate  images  of  the  target  locations.  Two  methods  of  image  generation  from  the  multi¬ 
static  data  were  developed  and  tested  in  computer  simulations;  (i)  image  formation  via  the 
classical  method  of  time-reversal  (or  field  back  propagation)  imaging  [2]  and  (ii)  use  of  the 
MUSIC  pseudo-spectrum  to  form  the  images  [4,  5, 6].  Both  of  these  methods  use  the  singular 
vectors  in  the  SVD  of  the  K  matrix  where,  however,  the  classical  time-reversal  scheme  uses 
those  vectors  associated  to  the  dominant  singular  values  while  the  MUSIC  method  employs 
those  singular  vectors  associated  to  the  small  singular  values.  It  was  shown  and  demon¬ 
strated  via  computer  simulation  that  the  classical  method  works  well  in  cases  of  a  few  well 
separated  targets  (well  resolved  targets)  but  fails  if  the  number  of  targets  is  large  or  they 
are  closely  located.  On  the  other-hand  the  MUSIC  method  works  independent  of  these  as- 
sumptions  and  returns  “super-resolution”  estimates  of  the  target  location  that  are  far  better 
than  would  be  obtained  via  the  classical  time-reversal  based  imaging  method. 

The  cmrent  final  project  report  included  a  discussion  of  the  effects  of  grating  lobes,  clutter 
and  additive  noise  on  the  image  formation  process.  It  was  argued  and  later  demonstrated 
in  the  computer  simulations  that  the  presence  of  a  inhomogeneous  medium  between  the 
phased  array  antenna  and  the  targets  can  significantly  reduce  the  grating  lobes  and  increase 
the  performance  of  the  algorithms.  Although  earlier  studies  [7]  have  shown  that  the  presence 
of  such  a  backgroimd  can  increase  the  resolution  of  the  time-reversal  imaging  process  the 
effect  demonstrated  here  appears  to  be  different  and  is  simply  due  to  the  randomizing  of  the 
grating  lobes  of  the  antenna  radiation  pattern  and  not  due  to  an  effective  increase  of  the 
aperture  of  the  antenna  array.  A  simulation  of  the  performance  of  the  MUSIC  algorithm 
in  the  presence  of  two  adjacent  and  dominant  targets  and  a  number  of  clutter  targets  and 
additive  noise  was  presented  and  illustrated  the  robustness  and  super-resolution  ability  of 
the  algorithm. 
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Both  the  classical  time-reversal  imaging  scheme  as  well  as  the  MUSIC  algorithm  require 
knowledge  of  the  backgrotmd  Green  function.  In  the  Tech  Sat  simulations  the  background 
was  modeled  as  a  thin  phase  screen  and  simulations  were  performed  for  cases  where  this 
background  was  known  and  also  for  cases  where  the  free  space  background  Green  functions 
were  employed  in  the  image  formation  process.  The  simulations  indicate  that  excellent  resiilts 
are  obtained  when  the  correct  Green  function  is  employed  in  the  image  formation  process 
but  that  performance  degraded  in  some  cases  significantly  when  the  back  propagation  or 
MUSIC  algorithms  were  implemented  using  the  firee  space  Green  function.  An  important 
avenue  for  research  is  the  estimation  of  the  background  from  the  multistatic  data  and  also 
firom  ancillary  data  obtained  using,  for  example,  the  phased  array  antenna  but  other  emission 
somrces  such  as  are  employed  in  Silent  Sentry  and  similar  projects.  Within  the  context  of  the 
Tedi  Sat  project  we  mention,  for  example,  that  the  reflections  off  the  ionosphere  which  are 
assmned  time-gated  out  in  the  simulations  presented  here  can  be  used  to  form  an  estimate  of 
the  index  of  refraction  distribution  of  this  medium  which  can  then  be  employed  to  generate 
an  estimated  background  Green  function.  An  inversion  algorithm  based  on  the  distorted 
wave  Born  approximation  has  been  recently  developed  [8]  for  precisely  such  cases  and  will 
be  integrated  into  the  time-reversal  based  schemes  discussed  here  in  the  futmre. 
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6  Table 


figure 

code 

OLk  in  m. 

Xm  in  m. 

Tm 

L 

1 

apOl.m 

-206,-104,-1,101,203 

- 

- 

20,  40 

2 

-206,-104,-1,101,203 

0 

1 

40 

3 

ap02.m 

-155,-104,50,75,153 

0 

1 

156 

4 

apOS.m 

-155,-104,50,75,153 

-201,-1,99 

1,1,1 

156 

5 

ap>04.m 

-155,-104,50,75,153 

-201,-1,99 

1,1,1 

156 

6 

ap04.m 

-155,-104,50,75,153 

-201,-1,99 

1,1,1 

156 

7 

ap04.m 

-155,-104,50,75,153 

-201,-1,99 

1,1,1 

156 

8 

ap04.m 

-155,-104,50,75,153 

-201,-126,-101,-96,-1,3,49,99 

.i,.i,i,i,.i,.i,.i,.i,.i,.i 

156 

Table  1:  Parameters  used  in  the  simulations 
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8  Figures 


Unstructured  Phased  Array 

•  Centimeter  wavelength 

•  5  to  15  independent  antennas 

•  kilometer  or  more  deployment  radius 

•  multistatic  data  acquired 

•  detection  and  tracking  of  MGT 

Clutter  suppression  performed  using 
time-gating  and  doppler  filtering 


Narrow  frequency  band 
Very  large  sparse  3D  array 
Intervening  ionosphere  (multipath) 


Figure  1:  Illustration  of  the  Tech  Sat  phased  array  system  for  MGT  detection  and  tracking 
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Figure  2:  Doppler  and  time-gating  are  used  to  reduce  the  size  of  the  ground  patch  and  thus 
reduce  clutter.  In  this  figure  an  exciting  voltage  ek{u})  is  applied  to  a  transmitting  antenna 
element  and  the  reflected  signal  from  the  target  generates  the  voltage  Vj{u})  at  a  receiving 
element.  This  received  voltage  is  input  to  a  Doppler  rejection  filter  to  generate  the  clutter 
reduced  signal  Vj{(j). 
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Transmitting  antenna 


Figure  3:  Within  the  Born  approximation  the  multistatic  data  matrix  K  can  be  expressed 
as  the  sum  of  outer  products  of  transmitter  (down  going)  gt  and  receiver  (up  going)  Qr  Green 
function  vectors.  The  Tm  are  the  effective  target  reflection  coefficients. 
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Computer  Simulation  Model 

Fresnel  Approximation 


smc[-^(®  —  rE„)] 


Thin  phase  screen  model 


,2tv 


T{x)  =  exp[i27rAo/(a;)];  f{x)  =  sin(— a:); 
Down-going  wave 

No. 

t/_(r,  r„)  =SxY^  T{xj)G{Tj ,  r„)G(r,  Tj) 
i=i 

Up-going  wave 

Na 

U+{r,  r„)  =SxYl  R{xj)U-{rj,  r„)(7(r,  r^) 

3=1 


Figure  4:  A  thin  phase  screen  model  is  used  to  simulate  the  effect  of  the  ionosphere.  This 
model  can  also  be  the  basis  for  a  parameterized  model  for  the  backgroimd  Green  function 
vectors.  In  the  figure  the  quantities  T  and  R  are  the  transmission  and  reflection  coefiicients 
of  the  thin  phase  screen  model  and,  as  a  first  approximation,  are  given  by  a  sinusoidal  model. 
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Figure  5:  (Top)  Plot  of  the  magnitude  of  the  eigenvalues  of  the  time-reversal  matrix.  (Bot¬ 
tom)  plots  of  the  (real)  phase  of  the  two  eigenvectors  corresponding  to  the  two  non-zero 
eigenvalues  (solid)  and  of  the  phase  of  the  comply  conjugate  of  the  two  Green  function 
vectors  (dashed). 
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Figure  6:  (Left)  Images  generated  by  the  two  eigenvectors  shown  in  Fig.  5  and  (right)  images 
generated  by  the  complex  conjugate  Green  function  vectors  shown  in  Fig.  5.  The  “X”  on 
the  images  indicates  the  location  of  the  targets. 
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Pseudospectrum  from  mus 


Figure  7:  The  pseudo-spectrum  computed  for  the  simulation  depicted  in  Fig.  6.  The  peak 
values  of  the  pseudo-spectrum  are  given  as  text  in  the  figure  and  indicate  that  exact  results 
were  obtained  for  both  the  x  and  z  location  estimates. 
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plot  of  absolute  value  of  eigenvalues 


phase  of  eigenvectors  and  Green  ftinctlon  vectors  across  array 


Figure  9:  (Top)  Plot  of  the  magnitude  of  the  eigenvalues  of  the  time-reversal  matrix.  (Bot¬ 
tom)  plots  of  the  (real)  phase  of  the  two  eigenvectors  corresponding  to  the  two  non-zero 
eigenvalues  (solid)  and  of  the  phase  of  the  complex  conjugate  of  the  two  Green  function 
vectors  (dashed). 
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Figure  10:  (Left)  Images  generated  by  the  two  eigenvectors  shown  in  Fig.  9  and  (right) 
images  generated  by  the  complex  conjugate  Green  ftmction  vectors  shown  in  Fig.  9.  The 
“X”  on  the  images  indicates  the  location  of  the  targets. 
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Pseudo-spectrum  from  mus 


Figure  11:  The  pseudo-spectrum  computed  for  the  simulation  depicted  in  Fig.  10.  The  peak 
values  of  the  pseudo-spectnnn  are  given  as  text  in  the  figure  and  indicate  that  exact  restilts 
were  obtained  for  both  the  x  and  z  location  estimates. 
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Figure  12:  Plots  of  the  field  intensity  over  the  phase  screen  (top),  over  the  ground  with  no 
phase  screen  present  (middle)  and  over  the  ground  with  a  phase  screen  present  (bottom). 
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Figure  13:  Plots  of  the  intensity  of  the  time-reversed  field  over  the  ground  with  no  phase 
screen  present,  with  phase  screen  present  and  with  phase  screen  data  but  free  space  back 
propagation. 
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Figure  14:  Same  as  for  fig.  13  but  with  a  different  antenna  geometry.  Note  the  larger  grating 
lobe  period  due  to  the  smaller  minimum  separation  between  adjacent  antenna  elements. 
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Figure  15:  Plots  of  the  intensity  of  the  time-reversed  field  over  the  ground  with  a  phase 
screen  present  for  the  cases  of  three  targets  located  at  Xm  =  -200,0,  and  -1-100  meters. 
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Plot  of  the  Pseudo-spectrum  with  no  additive  noise 
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Plot  of  the  Pseudo-spectrum  with  additive  noise 


100 
60 
60 
40 
20  ‘ 
0 


-0.25  -0.2  -0.16  -0.1  -0.05  0  0.06  0.1  0.15  0.2  0.25 


Figure  16:  Plots  of  the  pseudo-spectrum  corresponding  to  the  same  conditions  and  set  of 
parameters  used  for  the  conventional  time-reversal  images  presented  in  15  but  with  no  phase 
screen  present.  The  top  plot  is  the  pseudo-spectrum  without  additive  noise  while  the  bottom 
had  additive  noise  with  a  noise  ampUtude  =  0.1. 
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Plot  of  the  Pseudo-spectrum  with  no  additive  noise 
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Plot  of  the  Pseudo-spectrum  with  additive  noise 
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Figure  17:  Same  as  in  fig.  16  but  with  a  phase  screen  present. 
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Plot  of  the  Pseudo-epectrumHiSing  phase  screen  steering  vector 
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Plot  of  the  Pseudo-spectrum^using  free  space  steering  vector 
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Figure  18:  Plots  of  the  noise  free  pseudo-spectrum  corresponding  to  the  same  parameters 
used  in  fig.  17.  The  top  plot  is  the  pseudo-spectrum  computed  with  the  correct  steering 
vector  while  the  bottom  plot  uses  the  free  space  steering  vector. 
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Plot  of  the  Pseudo-spectrum  with  no  additive  noise 
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Plot  of  the  Pseudo-spectrum  with  additive  noise 
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Figure  19:  Plots  of  the  pseudo-spectrum  resulting  for  two  real  targets  and  eight  clutter 
targets.  The  top  plot  is  the  pseudo-spectrum  for  the  no  noise  case  while  the  bottom  plot  is 
the  pseudo-spectrum  for  an  additive  noise  with  amplitude  Ai  =0.1. 
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